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Abstract 

We investigate both numerically and analytically the convective instability 
of granular materials by two dimensional traffic equations. In the absence 
of vibrations traffic equations assume two distinctive classes of fixed bed so- 
lutions with either a spatially uniform or nonuniform density profile. The 
former one exists only when the function V{p) that monitors the relaxation 
of grains assumes a cut off at the closed packed density, p c , with V{p c ) = 0, 
while the latter one exists for any form of V. Since there is little difference 
between the uniform and nonuniform solution deep inside the bed, the con- 
vective instability of the bulk may be studied by focusing on the stability of 
the uniform solution. In the presence of vibrations, we find that the uniform 
solution bifurcates into a bouncing solution, which then undergoes a super- 
critical bifurcation to the convective instability. We determine the onset of 
convection as a function of control parameters and confirm this picture by 
solving the traffic equations numerically, which reveals bouncing solutions, 
two convective rolls, and four convective rolls. Further, convective patterns 
change as the aspect ratio changes: in a vertically long container, the rolls 
move toward the surface, and in a horizontally long container, the rolls move 
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toward the side walls. We compare these results with the those reported pre- 
viously with a different continuum model by Hayakawa,Yue and Hong[Phys. 
Rev. Lett. 75, 2328 (1995)]. Finally, we also present a derivation of the traffic 
equations from Enskoq equation. 

P.A.C.S. numbers: 47.20.-k,46.10.+z,81.35.+k 
I. Introduction 

This paper is concerned with the numerical as well as the analytical analysis of the two 
dimensional traffic equations, termed model B in the literature [1] with an aim to understand 
the convective instability of granular materials. It has been known since Faraday [2] that the 
granular materials in a vibrating bed develop permanent convective rolls when the strength 
of the vibration exceeds a critical value. Unlike the well known Rayleigh-Bernard thermal 
convection in fluids, however, the origin of the granular convection has remained relatively 
unexplored since its discovery, but recently two simultaneous push from experimental side 
[3,4] with the use of MRI or X-ray method as well as from computer simulations based 
on the distinct element method [5] have substantially aided our understanding through 
visualization. Yet, the theoretical efforts to uncover the basic mechanism of this granular 
convection have not been remarkable, still largely focused on producing convective patterns 
through computer models and simulations. While many impressive experimental data are 
currently being piled up, theoretical development [6-9] in this area seems to be still far from 
complete. In order to have a deeper understanding of the granular convection, we find it 
essential to come up with a reasonable continuum model that contains a minimum number 
of control parameters yet captures some of the essence of granular convection. The goal 
is certainly to explore analytically and numerically the instability mechanisms that lead to 
many fascinating complex nonlinear dynamics behaviours. 

Our aim here is two fold: we first present such a continuum model along with previously 
unpublished numerical results based on two dimensional traffic equations, and second, we 
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carry out the linear stability analysis of the traffic equations and uncover the mechanism of 
granular convection as a supercritical bifurcation of a bouncing solution. This mechanism 
for the granular convection is essentially identical to the previously obtained scenario for 
a different continuum model, termed model A [8]. However, detailed simulations of these 
two continuum models reveal distinctively different convection patterns for different geome- 
tries and parameter ranges. Most notable departure between the two is the migration of 
convective rolls toward the surface and toward the walls in model B, as the aspect ratio 
changes. This is in sharp contrast with the results of the model A that predicts a series of 
rolls in the container. Further, the simulation results of the model B yield rich dynamical 
behaviors as we probe deeper into the parameter space, yet the analytical structure of the 
nontrivial patterns seem to be quite complex and require further extensive studies in the 
future. Our analytical study here focuses exclusively on the stability analysis of a uniform 
bouncing solution, but there exists a second class of solutions that are spatially nonuniform. 
The stability analysis of the latter is highly nontrivial and we only present a brief description 
of it in the Appendix B. Fortunately, however, since there is little difference between the 
uniform and nonuniform solutions deep inside the bed, we expect that the uniform bouncing 
solution may capture the essence of the bulk instability. For details, see the text. 

This paper is organized as follows. We first define the model equations in II and make 
attempt to derive the traffic equations from the Enskoq equation in the Appendix A. We will 
then present numerical results in II, and provide some insight into the stability mechanism 
of the granular convection in III. We also discuss several unresolved questions in IV. We now 
turn to the main text 

II. Two Dimensional Traffic Equations and Granular Convection 

Apart from its wide application to traffic flow problems [12], the traffic equations [13-15] 
have recently been proposed, as an alternative continuum model to Navier-Stokes equation, 
in investigating a variety of dynamical responses of granular materials. Examples include the 
granular relaxation under repeated tapping [16], density waves and jamming and clogging 
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phenomena [11,14,17], 1/f spectrum of hopper flow [18], and granular segregation [19]. But 
we point out that, unlike the studies to be presented in this paper, almost all the studies 
thus far have focused exclusively on one dimensional traffic equations [13-19], which certainly 
will not be able to describe the granular convection. Hence, we first present here the two 
dimensional version of the traffic equations [1], termed model B in ref.[l]: 

dp dpv x dpv z _ Q 
dt dx dz 
dv x dv x dv x 2 dp d 2 v x d 2 v x 

-m +v ^ + v ^ = - c » ,p d-x + "< a*r + ^ (2) 

dv z dv z dv z V(p,t)-v z 2 dp d 2 v z d 2 v z 

-m +v ^ + v ^ = — ; cJp d- z + ^ + ^ (3) 

where c 2 ~ T e is the sound speed, and p is the viscosity. The difference between the model 
B [1] and the model A [8] is the presence of the relaxation term in the z direction, which 
is represented by an average function V(p, i) with the relaxation time r. The net effect is 
for the void(or particle) to adjust its speed, v z , around the average value V(p) in a given 
time r. The origin of such a term has been discussed in Hong et al [16] in an attempt to 
introduce correlations among grains or voids in the diffusing void model(DVM)[20]. In the 
DVM, the void speed is only a function of local void density, namely v z = V(p) + diffusion 
term. However, a void is a compressible hydrodynamic object that changes and adjusts its 
shape to conform to the surrounding, not instantaneously, but in a given finite time, r. So, it 
may be appropriate to write down the time dependent equation for the velocity in a manner 
given by (3) than simply assuming a fixed value at a given local density. The presence of 
such relaxation process may be effectively equivalent to assuming a drag force acting on a 
void. 

The coupled equations (1) -(3) are fairly similar to the two-phase model [26] of a fluidized 
bed that has been widely used for mixtures of gas and granular particles. As shown in the 
Appendix A, functions in the traffic equations (l)-(3) may be inferred from the Enskoq 
equation [21]; namely —v z /r is the drag term imposed externally on the particle, or by 
vibrations. In the case of no interstitial fluid, its origin lies in the frictions of the front and the 
rear glass of the container and from the wall. Further, the Enskoq pressure, Tp{l + f(p)p/2) 



with f(p) the correlation function, produces an extra term V(p) in addition to the hard 
sphere pressure Tp. In this case, we make an important observation that the coefficient of 
V(p) is proportional to the gravity g, which will enable us to incorporate the vibrations of 
the bed. This observation makes sense, because the strength of the mean speed (also termed 
the drift velocity) is determined by the gravity and thus it is quite physical to assume that 
the mean speed of the void is also proportional to the gravity as demonstrated in [16] and 
Appendix A. In the moving frame of reference of the vibrating bed, the mean speed V(p) 
depends not only on the density but also on the way the box is shaken. In the moving frame, 
the effective gravity, g = —g + Auj 2 cos(ujt), is the only time dependent factor involved in 
the mean speed V(p). Therefore, without the loss of generality, we expect that the mean 
speed of the granular flow along the z direction assumes the following form: 

V(p) = V (p) ■ g = Vi{p) • (-1 + Tsin{ut)) (4) 

where Y = Au 2 /g with A the amplitude of the vibration, and u the vibrational frequency. 
Vo(p) is the mean speed in the absence of the vibration, such as the case in the hopper flow. 
In one dimensional traffic equation, we have used [16]: 

V (p) = a/(l + (3p 2 ) (5) 

Note that the functional form of V Q is not unique, and other exponentially decaying functions 
have also been proposed [13]. The V (p) assumed above has a long tail when the density 
increases, which might be unphysical since there should be a critical density, p c , for any 
granular material, at which all particles will be locked in if the density exceeds that critical 
density. This critical density may be proportional or equal to the closed packed density. 
Thus, we have employed a different form in our numerical simulations that assume the form 
(Fig.la): 

V(p) = (p c - pf9(p c - p) (6) 

where we have used a step function 9{p c — p) = when p > p c . A physically reasonable 
choice of (5 may be such that V'(p) < oo. (For j3 < 1, V\p) diverges at p c , and thus in our 
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stability criterion, we set the density of the uniform solution slightly less than the closed 
packed density.) Before investigating eqs.(l)-(3) numerically, we first check the consistency. 

A. Consistency check: static solutions without vibrations. It is worthwhile to 
examine a few special cases not only for the sake of general limitation, but also to check 
the numerical accuracy. When T = Auj 2 / g = 0, which corresponds to the case without the 
vibration, which we also term "the fixed bed solution", we expect a static solution where 
v = in all momentum equations. We note that there are two classes of static solutions. 

First, a spatially nonuniform solution. Regardless of the functional form of V(p), we 
find by inspection that there exists a spatially nonuniform static solution that satisfies the 
following differential equation: 

V (p) = JVp (7) 

This represents a situation where particles pile up from the top to bottom with a certain 
density distribution. For simplicity, the mean speed is assumed to satisfy eq. || The mean 
speed as a function of density is plotted in Fig. la. When the density of the granular particles 
exceeds the critical density, say for example p c = 0.501, the particles become locked and have 
a steady-state density p = p c . If we solve the eq. [F] and plot the density profile as a function 
of z (Fig. lb), we can see that the uniform density profile is created close to the bottom 
of the plate, and then the density decreases linearly to zero as we move to the top of the 
container. Physically, this is acceptable if we notice that when the particles neat the bottom 
are locked, they act more like a solid. On the other hand, those particles near the surface 
are rather loosely packed and ready to move, which is more like a traditional fluid, where 
the density linearly decreases as the altitude increases. Note that the granular bed is not 
subjected to vibrations at this time and thus this is a fixed bed solution. Furthermore, the 
total number of particles is conserved by the mass conservation law of Eq.|I[ We emphasize 
that the non-slip fluid boundary conditions have been imposed in solving the equations 
numerically, namely, v ± = and v» = [8] 
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Fig.l. (a) Mean speed V (p) as a function of density is plotted above(a), with the cutoff 
density p c = 0.501 and (5 = 2/3. (b) Static solution is solved analytically and the density, 
p, is plotted as a function of the height z, where z = is the bottom of the container. The 
numerical result is consistent with this, to a very high accuracy. 

Second, a spatially uniform solution. There is, however, a second class of solution that 
exists only when V(p) has a cut off. This solution is a spatially uniform solution with 
the density given by the closed packed density, p c , everywhere with V(p c ) = 0. Numerical 
simulations of the traffic equations reveal that it is the spatially nonuniform solution that 
is realized in simulations. But the stability analysis of the uniform solution clearly shows 
that this solution is stable, too, and thus it must be realized with presumably suitable 
boundary conditions. Further, deep inside the bed, there is very little difference between 
the uniform and nonuniform solutions, and thus the convective instability inside the bulk 
may be investigated by focusing on the stability of the uniform bouncing solutions. In the 
absence of the rigorous stability analysis of the nonuniform solution, this is only our guess. 
Future work must focus on the rigorous stability analysis of the nonuniform solution(see 
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Appendix B) and compare its result to those of the uniform solution obtained in section III. 
Before presenting the stability analysis, we first present numerical results. 

B. Numerical results: We now investigate the time dependent phenomena numerically. 
Numerical solutions are obtained by two different algorithms for the rectangular boundaries. 
The first one is the center-difference method, which has the accuracy of the second-order in 
space and first order in time. The second method employed here is the RUNGE-KUTTA 
algorithm which has the accuracy of the second-order in space and third order in time. 
The results show that both methods give consistent conclusions for the same boundary 
conditions. For most of the cases, we set the grid number equal to N x = 32 and N z = 32 
in order to improve the accuracy while accommodating the CPU cost. Lower or higher grid 
numbers have been examined as well, but no significant changes have been found. Initially, 
the uniform solution is assumed, and typically four cycles are taken to calculate the average 
of the convection after the stable periodical solution appears at each site in the lattice. In 
the stream line plot, all of the vector flow arrows are normalized by the maximum speed. 
Numerically, we observe two distinctive classes of solutions: the first one is the bouncing 
solutions, and the second one the convective rolls. We first examine the bouncing solutions. 
Note that the set of parameters that correspond to physical situations might be uj/2it = 
20Hz, T e ~ 3, and ji — 1/R as 0.5 with the Reynolds number R — UL//j, s rs 2, because 
the kinetic viscosity of the granular material, /i s ~ 5 • lCT 3 m 2 /s and the typical velocity 
U ~ 10cm/ s [8]. For pure numerical reasons, however, in what follows, unless mentioned 
otherwise, all the simulations are conducted with the parameters, /i — l,cu — 1,t — 1 and 
T e = 2, L=32 with V Q (p) given by Fig.l. 

(a). Bouncing solutions with T < T c : For a fixed set of control parameters, bounc- 
ing solutions, first termed in [8], appear when the vibration strength T is less than the 
numerically determine critical T c 1.5. Upon increasing T, this bouncing solution becomes 
unstable and bifurcates into a different branch and convective rolls appear. We have exam- 
ined the case where the vibrational acceleration is less than the gravity, in other words, the 
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dimensionless variable T is less than one, which is certainly the regime below the critical 
strength. In real situation, for T < 1, nothing much really happens inside the box except 
perhaps the volume decreases due to settling caused by vibration. However, in our simula- 
tion, due to the presence of the top layer, we expect the bouncing solution to appear even for 
T < 1. We have performed the simulation using different initial conditions such as a uniform 
density distribution, a linear density distribution, or a fluctuation in the velocity field. For 
all cases, we have found that after a relatively long CPU time, the convection pattern disap- 
pears, which implies that the convection does not persist in this parameter regime. However, 
if we measure the local density or velocity at different locations, we clearly observe oscilla- 
tory solutions that persist everywhere inside the container. This is termed the "bouncing 
solution" that exists before the onset of the convection: the bulk of the granular material 
oscillates up and down, like a single solid ball, with a central frequency of the vibrating 
container and some harmonics. The harmonic oscillations of the density and velocity at the 
same lattice site are displayed in Fig. 2a and 2b. 
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Fig. 2. Bouncing solutions for (a) the velocity v z and (b) the density at the same fixed 
location (x,z) = (15,15), where no convection rolls appear and the bulk of the particles moves 



as a single particle. T < T c . 



(b). Convective domain with T > T c : When T > T c , the bouncing solutions disap- 
pear and the permanent convective rolls appear inside the bulk. In order to view this slow 
convection, we have to take the average of the velocity field over many periods. Typically 
four cycles are used, since the amplitude of the convection is relatively small in comparison 
to the amplitude of the vibration. This is true in experimental systems as well. In the 
stream-line plot, we have normalized the velocity field, by taking the maximum velocity 
component as the maximum length of the arrows. 

In Fig. 3a convection patterns are illustrated for T = 2.25 > T c by plotting the streamline 
of the velocity field, where V = (p c — p) /3 0(p c — p) with f3 = 2/3 and p c = 0.501 as given in 
Fig.l, is used. Notice that the flow direction close to the wall is toward the upward direction, 
and at the center of the container, the particles try to move up. The density profile is also 
plotted in Fig.3b, where the average of the density over one period remains close to the 
static solution. At any given time the density profile varies around such an average velocity. 
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Fig. 3. Convection pattern and density profile: (a) When T > T c , convection patterns 
are illustrated by plotting the streamline of the velocity field, where the velocity field is 
normalized by the maximum velocity component v z , with Vq given by Fig.l. (b) The density 
profile as a function of position. The average of the density over one period remains close to 
the static solution(+) and at any given time the density profile (circle) swings around such 
an average velocity. 

Next, it has been found that the geometry plays an important role in the convection pat- 
tern formation [7,22] and we examine its effect on the convective instabilities. For example, 
the change of the aspect ratio(defined as width/height) may affect the relative positiona, 
the number of convection rolls, and the intensity of the convection with respect to the same 
vibration parameters. We have examined two particluar cases: 

First, the numerical results displayed in Fig.4a and 4b show that the location of the 
convection rolls migrate when the aspect ratio changes. For the aspect ratio, Nx/Nz = 
0.5, convection rolls move toward the top of the container, and the particles close to the 
bottom are almost locked allowing no relative motion to exist. This is consistent with 
the experiments and MD simulation results. Convection rolls are shown in Fig.4a and the 
corresponding density profile is plotted in Fig.4b. Note that the stability analysis of model 
A [8] predicts a series of rolls when the vertical axis increases as a multiple of the horizontal 
width. The form of V(p) with a cut off at the closed packed density seems to be crucial in 
capturing the real convective patterns. Further, the average density curve (+) is very close 
to the static solution without the vibrations. The second curve (— ) is the density profile 
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at the stop time t, which clearly shows that the density profile varies in time around the 
equilibrium density profile. 
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Fig. 4. Convection pattern for a small aspect ratio N x /N z = 0.5. All other parameters 
remain the same as in the case of N x /N z = 1.0. (a) The two convection rolls move to the 
top region of the container and most of the particles closed to the bottom are locked during 
the whole vibration cycles. This can be verified in the density profile at a fixed point x=16 
(b), where the time dependent density curve(circle) varies around the static solution(cross), 
but much less at the bottom than at the top regime. 



Second, we have also examined the case with the large aspect ratio, N x /N z = 2.0, (results 
are shown in Fig. 5a), where the width is twice larger than the height. Here, four convection 
rolls appear with two small ones near the top corners (Fig. 5a). If we extend the width even 
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longer to the aspect ratio, N x /N z = 4, the convection pattern will only exist around the 
narrow regime of side-boundaries, represented by the zero intensity of convection (Fig.5b). 
It seems that convection rolls only appear close to the boundaries, which indicates that the 
friction on the wall presumably plays a significant role in the origin of the vibration-induced 
convections. Such observations are consistent with the observations made by Taguchi in 
Molecular Dynamics simulations [22]. 
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Fig. 5. Convection pattern for large geometry ratio, (a) Convection pattern is drawn for 

a large aspect tratio, N x /N z = 2.0. All other parameters remains the same as in the case 
of N x /N z = 1.0. Two main convection rolls appear at the center, with another two visible 
small convection rolls at the top corners of the container. The two small convection rolls 
have opposite spins relative to the major ones, (b) Two convective rolls for the aspect ratio 
N x /N z = 4.0 are illustrated for the same parameter regimes. The two rolls exist near the 
walls. They are separated by a wide non-convection regime, where the velocity field equals 
zero. 
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C. Pattern dependence on V (p) While we have not systematically examined the stability 
diagram of the many dimensional parameter space, we report here a few cases simply to show 
the richness of the problem associated with the traffic equations. In the previous arguments, 
we have only assumed that Vo(p) likely has a cutoff when the density exceeds a critical 
density (at which the particles become locked), and reaches the maximum when the density 
approaches zero (since there are no interactions). On the other hand, it is also reasonable 
that Vo(p) might have a long tail in the high density regime, which may represent very rough 
grains or grains non-uniform in size. Hence, we have examined numerically the case with 
the following V (p). 

V (p)=a/(l + f3p 2 ) 

Our numerical results indicate that the direction of the convection rolls is quite different from 
the results obtained with a V (p) that assumes a cutoff. Moreover, the four roll convection 
pattern appears when the long-tail Vo(p) distribution is assumed. In Fig. 6 is plotted the 
convection streamline with a = 0.2, (5 = 0.5. As we increase the vibrational frequency, the 
high intensity regime of the convection rolls moves toward the boundary. At some parameter 
regime, the two roll patterns evolve to four rolls. In Fig. 7, the four roll convection pattern 
is plotted. 

In Fig. 6, Fig. 7, Fig. 8, and Fig.9, a wide parameter range has been examined. It seems 
that the appearance of the two roll pattern or four roll pattern not only depends on the 
aspect ratios, but also on the vibrational frequency and amplitude. Future work must 
focus on exploring the parameter space as well as the effect of boundary conditions to the 
convective patterns. We now present the linear stability analysis of the two dimensional 
traffic equations. 
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Fig. 6. Convection patterns for slow mean speed: For slow mean speed a = 0.1, vibrating 
amplitude A = 2.0, (3 = 0.5. Two convection rolls appears for the different vibrating 
frequency, where higher frequency corresponds to a higher flow rate close to the boundaries, 
respectively from left to right uj = 1.0 (a) and uj = 15.0 (b). 
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Fig. 7 Convection pattern with four rolls: (a) Four convection roll appears when we 
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increase the amplitude of the mean speed a = 0.5 with the frequency cu = 1.0. (b) As we 
increase the vibrating frequency, two roll convection patterns appears with a higher flow 
rate close to the boundaries. 
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Fig. 8 Convection pattern with four rolls: (a) Four convection roll persists when we 
increase the amplitude of the mean speed a = 5.1 with the frequency uo = 1.0. (b) As we 
increase the vibrating frequency, two roll convection pattern appears with a higher flow rate 
close to the boundaries. 
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Fig. 9. Convection patterns with high frequency: (a) For the same Y = 2.0 as in Fig.8, 
which corresponds to the induced vibrating gravity acceleration, (b) Four rolls convection 
pattern disappears when we increase the frequency to uj = 15.0. Instead, the time average 
of the one period of oscillation is very random and the magnitude of the convection decays 
as the function of the periods, where in the calculation, a = 0.5, (5 = 0.5. 



III. Linear Stability Analysis of Traffic Model 

Before studying the time dependent properties of the vibrating beds, any model must 
make sure that the fixed bed solution, which is the ground state at T=0 or equivalently 
r = [23], must be stable. 

(a) Fixed bed solution: Eqs.(l)-(3) assume two different fixed bed solutions. The first 
one represents a uniform granular block characterized by the homogeneous density, p = p c 
such that V(p c ) = 0, and v = 0. But there is a second non-uniform solution, given by eq.(7). 
Both solutions represent fixed bed solutions. We only consider the stability of a uniform 
solution. Let p = p c + Pi and v = v x . Substituting these into eqs.(l)-(3), we find that p x 
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satisfies the second order equation in time: 



Pi + Pi/t - pV 2 pi - T e V 2 pi = 



(8) 



Assuming, p\ = p q (t)sin( / rrmx)e' iqz , we find: 



p q + B (q)p q + D (q)T e p q = 



(9) 



with 



B (q) = 1/r + p(n 2 m 2 + q 2 ) 



(10a) 



D (q) = Te^ 2 m 2 + q 2 ) 



(106) 



For p q (t) ~ e CTt , it is straightforward to show that Rea is always negative, i.e,Rea < 0. 
Thus, the uniform fixed bed solution is stable. For the stability analysis of a nonuniform 
solution, see Appendix B. 

(b) Stability of the bouncing solution: When the system is subjected to vibrations, how 
do uniform and nonuniform solutions evolve to eventually produce convective rolls? The 
standard way of proceeding the stability analysis is first to find the basic solution, and 
then exmaine its stability by performing the linear stability analysis of this basic solutions. 
The onset of convection may be identified as a point at which such basic solution becomes 
unstable. In the thermal convection problem such as the Rayleigh- Bernard convection, 
the basic solution is the laminar solution with a temperature dependent density profile 
[25]. Since there exist two fixed bed solutions, we need to study the stability of these two 
solutions separately. In what follows, we present only the stability analysis associated with 
the uniform fixed bed solution, and show in the Appedix B why the similar analysis for 
the nonuniform solution is nontrivial. Our analysis in what follows is based on the belief 
that since there is no difference between the two deep inside the bed, the mechanism for the 
convection for the uniform and nonuniform solutions might be the same. When r ^ 0, the 
basic solution associated with the uniform fixed bed solution is also a homgeneous solution in 
density, i.e. p = p Q , which then undergoes harmonic oscillations. This represents a uniform 
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granular block, represented by a single ball [25], moving up and down bouncing at the plate. 
We previously termed such a solution, a bouncing solution and such a bouncing solution can 
be readily found by inspection of the traffic equations (l)-(3). Its horizontal velocity is zero, 
but the vertical velocity, u (t), is given by 

v z = uj (t) = —V {p ) + V{p , r, r)(sin(ut) - ujTcos(ujt)) (11) 

where V(p ,T,T) = TV (p )/(l + uo 2 t 2 ). The numerically observed bouncing solutions in 
Fig. 2 may represent such a uniform bouncing solution. 

We now examine the bouncing solution in some detail, which is essentially identical 
to a single ball dynamics, if we focus on the motion of the center of mass of the block. 
The dynamics of a ball on a vibrating plate is non-trivial for high vibration strength, but its 
trajectory obeys the simple Newtonian dynamics for low Y with basically two characteristics: 
(i) the ball moves together with the plate for some time, (ii) the ball is launched to the space 
at a particular time t that depends on the vibration strength. We assume that at the time 
of launching, the relative velocity between the block and the plate is such that the relative 
distance, A(t), between the ball and the plate starts to increase at the time of launching. 
Hence, t a is determined by the condition: 

r 

—(sin(ujt ) - uTcos(ut )) = 1 + e (12) 

1 + uo z r z 

where e is determined by imposing the self-consistency A(t = t +) > (eq.14). Note that 
(12) ensures that 

(1 + e + cuV) 1 / 2 < T (13) 

For T to be order one, r needs to be of order ~ 1/ou. Next, since the velocity is given, we 
find that the position of the ball in the moving frame, A(t) = J t * uj(t')dt', is given by: 

TV 

A(t) = -V {p ){t - t ) + —((cos(ut ) - cos{ut))/u - r(sin(ut) - sin{u:t )) (14) 

1 + LO Z T Z 

In Fig. 10 are shown typical v z (t) and A(t) as a function of time. Note that A is the 
relative distance between the ball and the plate and thus cannot be negative: A(t) is positive 
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for t Q < t < ti where 0J o (t o ) = A(t ) = and A(ti) =0. Tf — t± — t a , which is a function 
of the vibration strength T, is the flying time of the ball. For a set of parameters used in 
simulations, i.e \i = uj = r = 1.0, and V (p ) = 0.5, we find (t a , ti) changes from (2.356, 4.867) 
with Tf = 2.51 for T = 2 to (1.197, 5.0368) with T f = 3.84 for T = 5. For larger T, the flying 
time Tf obviously increases. 

1.5 

1 

0.5 


-0.5 

-1 
-1.5 

-2 

1 2 3 4 5 6 7 8 9 10 11 12 

Fig. 10 Bouncing solutions. The velocity,u; (t) (lower one), and the position(upper one) 
of the ball relative to the plate as a function of time. The ball is launched at t — t D (Eq.l4) 
with a positive velocity v z (t ) > 0, and flies in the air until it hits the plate again at t — ti 
with A(ii) =0. The ball stays at the plate and is relaunched at t = t a + 2-n/uo. 

Let us now examine the stability of this bouncing solution. Let p = p + Pl and v = 
u) {t)z + vl and substitute p and v into equations (l)-(3). Taking only the linear term, we 
find: 

dtpi + PoV ? • v L = (15a) 
d t v L = -V^p L - — + z^^-p L (-l + Tstnicut)) + pV 2 v L (156) 

Po T T 

where £ = z- fo^(t')dt' = z- A(t). £ measures the point of the granular block in the box 
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frame that undergoes vibrations. Now, taking the divergence of eq.(15b) and substituting 
p Q Vg • v l = —d t pL, we find the equation for p L : 

(L + L x )p L = (16) 

where 

L = d 2 - T e V\ - iiVp t (17) 

L x = r-% + p T- x V' {p ){-\ + rsm(wt))^ (18) 

In order to study the stability of the bouncing solution, we now seek a solution of the 
form: 

Pl = p L , m , q (t)^mx/L)e^ +A ^ (19) 
Note that m is integer and £ + A = z. Substituting (19) into (16), we obtain: 

p q + B{q)p q + iC(q)p q + D(q)p q + iE(q)p q = iL q (t)p q (t) + M q (t)p q (t) + iN q (t)p q (20) 

where 

B(q) = B (q) = p(fr 2 m 2 + q 2 ) + 1/r (21a) 

C(q) = -2qV (p ) (216) 

D(q) = D {q) - q 2 (V 2 ( Po ) + ^(1 + u 2 r 2 ) (21c) 

£(g) = qp(7r 2 + q 2 ) - (21<Q 

r 

L ? (t) = -2gn>(sm(o;t) - cjrcos(cut)) (21e) 

2 2 i 

M,(i) = g 2 [ - T ~ T 2 V 2 cos(2cut) - urT 2 V 2 sin{2ut)} (21/) 

N q (t) = -T^^stn(ujt) (21g) 

where % = n/L with L the box size and p q (t) = pL,m,q(t). For the most unstable mode, we 
set m=l. We now change the time, t t Q + T and focus on the one oscillation from T = 
to T = 2ti/uj. Then, (21e)-(21g) are replaced by: 

L q (t) — > — 2qTV[{sin{ujt ) — uoTCOs{ujt ))cos{uoT) + (cos(ut ) + UTsin(Ljt ))sin(uT)] 
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M q (t) -> g 2 r 2 T> 2 [( - T - cos(2a>t ) - wrsin(2wt ))cos(2wT) 

2 

CO^T 2 1 

— ( sin(2cut ) + uTcos(2ut ))sin(2u>T)] 



N (t) -> - Fp ° gV ° (sin{ut )cos{u}T) + cos{u}t )sin(u}T)) 

T 

The stability of eq.(20), is highly nontrivial because it contains the time dependent 
inhomogeneous term, and we employ here an approximation method that was first used in 
ref.[8], namely we replace the time dependent function L q (T), M q (T) and N q (T) by their 
average over one period, namely < L q (t) >= TJ 1 J^ f dtL q (T) with Tf = (ti — t Q ), where 
the time dependent parts in L q (t), M q {t) and N q {t) in Eq.(21) are replaced by their average 
values over the flying time: 

sin(nuT) — > — ^r(l — cos(ncuTf)) 

cos(nuT) — > — \—sin(nujTf) 

where n=l or 2. The validity of such an approximation may be justified since we are 
interested in the long time behavior and the system continuously oscillates in time. Hence, 
we have replaced the inhomogeneous eq.(20) into the homogeneous one: 

p q (T) + B(q)p q (T) + i<3(q)p q (T) + D(q)p q (T) + iE q {T)p q {T) = (22) 

where 

C(q) = C(q)- < L q (T) > (23a) 

D(q) = D(q)- < M q (T) > (236) 

E(q) = E(q)- < N q (T) > (23c) 

The stability of the homogeneous equation (Eq.22) with time independent coefficients 
is easily determined by a standard linear stability analysis, namely we set p q (t) = e at and 
substitute this to (22) to obtain: 

a 2 + {B{q) + iC{q))a + D(q) + iE{q) = (24) 
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The condition for the instability is when the perturbations grow in time, i.e., Rea > 0. Note 
that the solution of Eq. (24) is given by 

°± = - m \ lC{q) ± \y/(B(q) + %C{q)Y - AD(q) - AtE(q) (25) 

In order to extract the real part, we notice that the argument of the inside the square 
root may be written as Re** with R = ((B 2 - C 2 - AD) 2 ) + ABC - 2E) 2 ) 1 / 2 and tari(p = 
2{BC - 2E)/(B 2 - C 2 - AD). The real part of a±: 



Re[*±\ = -f ± \\[^ - Df + (E- \BCf + - D\ (26) 

When the real part of eq.(26) becomes positive, the bouncing solution becomes unstable. 
Hence, the instability condition for the most unstable branch becomes, Re(a + ) = —B/2 + 
R l l 2 cos(<j)/2) > 0, which is rewritten by the following inequality: 

(BC-2E) 2 > B 2 (B 2 - \B 2 -C 2 - AD\) (27) 

Since the granular block is only in motion for the flying time, < T < Tj, but sits on 
the plate for Tf < T < 2n/uj, the stability of the granular block needs to be determined by 
the effective instability condition, namely for the flying time, Eq.(27) is operative, but for 
the rest of the time, the coefficients for the fixed bed solutions, (10a) and (10b), must be 
used. We thus define the effective instability condition: 

a eff = T f [(BC - 2E) 2 - B 2 (B 2 - \B 2 - C 2 - AD\) 

+ {2-k/uj - T f )(-B 2 (B 2 - \B 2 - AD Q \) > (28) 

When a e ff > 0, the bouncing solution becomes unstable. Since B(0) = fin 2 + 1/r, 
(7(0) = C(0)- < L q (0) >= 0, D(0) = D o (0)- < M q (0) >= D o (0) = T e rc 2 , E(0) = 
E(0)- < N q (0) >= 0, we find a eff (q = 0) = -^(O) 2 ^^) 2 - \B 2 (0) - AD o (0)\)f . If 
B 2 (0) > AD o (0) or (^tt 2 + 1/r) 2 > 4T e vr 2 , then a eff (0) < 0. If B 2 (0) < AD o (0), then 
fJ e //(0) = 2 J B 2 (0)( J B 2 (0)-2D o (0)) can be either positive or negative. If a eff (0) > 0, then the 
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bouncing solution is always unstable. This determines the criterion for the existence of the 
bouncing solutions in the traffic equations: (/i7r 2 + 1/r) 2 < 4T e 7r 2 and (pit 2 + 1/r) 2 > 2T e TX 2 . 
And thus, if 

(/i + 7r 2 + l/r) 2 <2T e 7T 2 

then, uniform bouncing solutions do not exist. Otherwise, we expect the bouncing solutions 
to appear. 

In Fig. 11 is displayed a e ff as a function of the wave number q for r = u,p — 1, p — 0.53, 
V (p ) = 0.5, V^(p ) = —10., L = 32, and T e = 2.0. Note that for T < T c , <7 e ffc is always 
negative, and thus the bouncing solutions are stable. As we increase the control parameter 
T, the peak of <7 e // moves up and becomes zero at the critical T c = 2.27, which signals the 
instability of the bouncing solutions, and thus is identified as the onset of convection. For 
T > T c , there is a band of wavenumbers where a e ff > 0. The maximum growth rate is 
q c — 1.11 for T = 2.32 and q c = 0.95 at T = 2.27, which are of order one. The corresponding 
wavenumbers for the convective rolls are: A = 2n/k ~ 6. Since there are two rolls in a box 
of size L=32, the wavelength of the convective rolls is about 8, and our analysis is consistent 
with the numerical results. For other parameters, the qualitative feature is the same, namely 
for T < T c , a e ff <) for all q, and at a critical T c , the local maximum of a e ff crosses zero at 
q = q c , and for T > T c , there is a band of wave numbers for which a e ff is positive. 
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Fig. 11 The effective instability condition cr e ff(q) (vertical axis) is plotted against the wave 
number q for different vibration strength T. At the critical r c = 2.27, the local maximum of 
the growth rate crosses the zero. For T = 2.2 < r c , a e ff < for all q, while for T c < T = 2.32, 
there is a band of wave numbers for which a e ff > 0. 

IV. Discussions 

We conclude this paper by briefly reviewing the present status of the convective instabiity 
of vibrating beds and simulations of the models presented here (Model B) and elsewhere [1]. 
From our analytical as well as numerical analysis, we confirm that the bouncing motion of 
particle collections plays an important role in both model A and model B. We have, however, 
several issues which should be resolved in the near future. 

First, the role of boundary conditions. We have employed no slip boundary conditions at 
the walls and the plate. Further, we have put the rigid walls at the top and thus suppressed 
the surface motion. Unlike the fluids, Granular materials have been shown to exhibit very 
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different motion near the wall in a zet flow experiment under gravity by Caram and Hong [20] 
and there has been some attempt to use the negative or positive slip to control the convective 
patterns [7] . More detailed studies to derive reasonable boundary conditions at the wall are 
required. 

Second, the role of interstitial fluid. The model A [8] assumes no interstitial fluid such 
as air, and it predicts a series of rolls for the vertically large aspect ratio. Model B on the 
other hand predicts that the convection suppresses in the bulk region but is confined near 
the surface, which appears to be in accordance with experiments and with the results of 
the two- fluid model with an interstitial fluid [26]. Perhaps, the origin of drag, whether it 
is coming from friction at the walls, or from the viscous effect of the interstitial fluid, may 
not be relevant once it is present. The suppression of the convection in the bulk is due to 
the locking mechanism of grains for near the closed packed density, which was taken into 
account in the present model B by a cut off in V(p), namely V(p) = for p > p c . The 
results of model B, in particular the migration of convective rolls toward the surface(Fig.4) 
and to the side walls(Fig.5) seem to be more realistic and closer to the experimental and 
MD simulation results than those of model A. We certainly need more extensive studies of 
both model A and model B to make quantitative comparison with experiments. 

Third, future work must focus on completing the stability analysis for the spatially 
nonuniform solutions(Appendix B) and its instability mechanism along with the nonlinear 
analysis. In addition, the study of the surface instability and its connection to the surface 
fluidization and excitations as well as the study of horizontal vibrations and the onset of 
liquefaction [9,27] remain important open problems, which may shed insight into the oscillons 
[28] and other related surface instabilities [29]. Such studies are currently under way and 
will be reported in future communications. 
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Appendix A: Derivation of Traffic eq. from Enskoq eq. 

We present here a brief derivation of the traffic equations from the Enskoq equation. We 
first recall how one obtains the hydro dynamic eq. from Boltzmann eq.[21] 



df 

-L + vV r f + F.V v f = J B (A-l) 

where F is the particle distribution function, and Jb is the collision integral. If x is summa- 
tion invariant, i.e., then xi + X2 = Xi + X21 where ' refers to the velocities of the two colliding 
particles after the collision, and / d 3 px(r, p) J B = 0. Now, one can derive hydrodynamic eqs. 
by multiplying / d 3 px(r,p) to (A-l) as detailed in Huang [30]. We now just focus on the 
Navier Stokes eq. 

(I- + u • V)u = F/m - -V(P - • u) + ^V 2 u (A - 2) 

at p 3 p 

For granular materials, we need to take into account two facts: first, particles are not 
point particles, but hard sphere particles with finite diameter D, so one should use Enskoq.eq. 
instead of Boltzmann eq [21]. Second, we need to take into account gravity and friction, for 
which the external force term F in (A-2) may be replaced by 

F/m = -gz - 711 (A - 3) 

The second term in (A-3) is due to the friction. Further, for Enskoq eq., there is collisional 
transfer. So, the pressure tensor has two components: 

p.. = pkinetic + pCT ^ _ ^ 

The lowest order to the collisional transfer is: P zz = x(p)p 2 /2 with x = (1 — p/2)/(l — p) 3 , 
where x is the two particle equilibrium correlation function at the contact point of two 
colliding particles. It is known that the Enskoq pressure has the following two components: 

P = T P (1 + Xp/2) = P ld eal + Pi (A- 5) 
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Let us focus only on z component. Then, the Navier-Stokes equation becomes: 

. d d . 1 d n . . . „. 

(— + -u 2 — )-u z = — (? — 7M Z — -r + viscous terms (A — 6) 

at oz pdz 

Now, 

dP/dz = T\p' + ^p 2 /2 + pp' X ] 
And the right hand side of (A-6) becomes: 

M-S = ~9-~- TWp/2 + p'x) - 
p dz 

Tdp | (V(p)-M z ) M-7) 



where r = 7 1 and 



\/(p) = --[/3 + x'p/2 + p'x] 

7 



Hence, we have shown that the z component of the Navier Stokes eq. reduces to the 
traffic eq. 

,d d , Tdp (V(p)-u z ) „ 2 , A ns 

(i* +u *ir) u * = — / + ^^ - + pV 2 u z (A -8) 

at oz p dz r 

For x and y components, simply use the Boltzmann eq. 

Note that within this derivation, the mean speed has the following form: 

g A-17p + 8p 2 
V{P) = -7 4(1 -pf 

with the coefficient proportional to the gravity. This is the reason for the functional form 
taken in eq.(4) in the text. For the equation of void, we need to change p — > (1 — p), for 
which V(p) decreases as a function of p. 



Appendix II: Stability analysis of a spatially nonuniform fixed bed solution 

We show here why the stability analysis of a spatially nonuniform fixed bed solution 
is complicated. In the absence of vibrations, T = 0, and the solution of this class, p Q (z), 
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is given by eq.(7). Since the solution is spatially inhomogeneous, for small T, we seek for 
bouncing solutions that are only a function of z: 

p{z)= Po {z)+p\z)s{z)e^ (B-l) 

v z (z) = v(z)e^ (B - 2) 

where p'(z) = dp(z)/dz. We put (B-l) and (B-2) into the traffic equations (l)-(3) and 
linearize to obtain: 

1 dp (z) . d 1 dp , . 

-ius(z) + (— H —)v(z) = (B — 3) 



p (z) dz dz p dz 

~ c2 °7 ^lT " iuv{z) " (1/r " ^ )v{z) = \ v (p°W {B ~ 4) 

Let Q(z) = ~j- o ^r = ~^V(p ). Then, we find eqs. for Q(z) and v(z): 

-iwQ(z)s(z) + (^--Q(z))v(z)=0 (B-5) 
dz 

We need to find the solutions for (B-5) and (B-6) that both satisfy the boundary conditions, 
s(z),v(z) — > as z — > oo, which appears to be nontrivial if we relax the condition that T is 
not a small parameter. Even though we are lucky to find the solutions, the stability analysis 
around these solutions still present difficulties. 
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